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A flexible spatio-temporal model is implemented to analyse ex¬ 
treme extra-tropical cyclones objectively identified over the Atlantic 
and Europe in 6-hourly re-analyses from 1979-2009. Spatial varia¬ 
tion in the extremal properties of the cyclones is captured using a 
150 cell spatial regularisation, latitude as a covariate, and spatial 
random effects. The North Atlantic Oscillation (NAO) is also used 
as a covariate and is found to have a significant effect on intensifying 
extremal storm behaviour, especially over Northern Europe and the 
Iberian peninsula. Estimates of lower bounds on minimum sea-level 
pressure are typically 10-50 hPa below the minimum values observed 
for historical storms with largest differences occurring when the NAO 
index is positive. 


1. Introduction. Extreme North Atlantic and European extra-tropical 
cyclones are a major source of risk for society. These natural hazards cause 
much damage and insurance loss in Europe due to extreme wind speeds/ 
flooding. Recent examples include the December 1999 windstorms Anatol 
and Lothar [Ulbrich et al. (2001)], and windstorm Kyrill in 2007 which 
resulted in large losses across most of Europe. Important scientific questions 
are as follows: 

1. How extreme (intense) can extra-tropical cyclones become? Or, more 
precisely, how much more extreme compared to the most extreme values 
recorded in short series of historical observations/analyses? 

2. How does the extreme behaviour vary spatially? 

3. How does the extreme behaviour vary in time due to modulation by 
large-scale climate patterns? 
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We consider sea-level pressure (i.e., cyclone depth) as a measure of cyclone 
intensity. Unfortunately, there are no simple physical arguments for how 
deep an extra-tropical cyclone can become. The most extreme events often 
deepen explosively with rapid decreases in central pressure, for example, 
storms known as bombs having pressure drops of more than 24 hPa in 24 
hours at 60 N. Explosive cyclogenesis depends on many factors, for example, 
the deepest recorded 20th century low of 913 hPa (the Braer cyclone of 
January 1993) deepened 78 hPa in 24 hours due to a combination of several 
factors such as available moisture and stratospheric conditions [Odell et al. 
(2013)]. The unlikely possibility that such conditions could be maintained for 
2 days gives a minimum value of SLP of around 990 — 156 = 834 hPa starting 
from a typical background state of 990 hPa. It should also be noted that 
SLP less than 650 hPa would correspond to mid-latitude geostrophic wind 
speeds faster than the speed of sound, which due to shock wave dissipation 
would be impossible to maintain energetically. In the absence of any more 
rigorous physical bounds, it is of interest to estimate bounds empirically 
using statistical approaches such as extreme value theory. 

Modelling cyclones poses an interesting challenge: the events occur ir¬ 
regularly in space and time with rates and magnitudes that are spatially 
heterogeneous and nonstationary in time (due to modulation by large-scale 
climate conditions). Furthermore, at any one location, very few extreme 
events are observed in short historical data sets. Here we model extreme 
North Atlantic cyclones using an extended version of the spatial point pro¬ 
cess model for extremes from Cooley and Sain (2010). The extension involves 
the inclusion of temporal covariates, the adaptation to irregularly occurring 
(i.e., random occurrence rather than fixed locations) extremes in space and 
the application to extra-tropical cyclones. 

2. Background and data. 

2.1. Extreme extra-tropical cyclones. There has been surprisingly little 
use of extreme value theory to investigate extreme cyclones [see Katz (2010) 
for a discussion about the lack of extreme value theory in climate science]. 
Lionello, Boldrin and Giorgi (2008) investigated changes in future cyclone 
climatology over Europe using the Generalised Extreme Value (GEV) dis¬ 
tribution to model pressure minima. Return levels were calculated over the 
whole North Atlantic domain without explicit characterisation of spatial 
or temporal heterogeneity. Della-Marta and Pinto (2009) and Della-Marta 
et al. (2009) used a Generalised Pareto Distribution (GPD) model to anal¬ 
yse future changes in extreme wind intensity. Three large nonoverlapping 
areas were considered, however, there was no formal consideration of spatial 
or temporal variation in the extremes. Sienz et al. (2010) used GPD models 
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extending the work by Della-Marta and Pinto (2009) to include temporal co¬ 
variates such as the North Atlantic Oscillation (NAO) and a linear trend but 
did not account for spatial variability. Bonazzi et al. (2012) used bivariate 
extreme value copulas to model the spatial dependence in footprints of peak 
gust wind speeds from a set of 135 damaging European cyclones. However, 
this study did not explicitly model the magnitude of many cyclones and so 
does not answer the question about upper bounds on cyclone magnitudes. 

2.2. Brief review of spatial extreme models. Davison, Padoan and Ri- 
batet (2012) identified three main classes of statistical models for spatial 
extremes: Bayesian hierarchical models (BHM), copula based models and 
max-stable process models. Although max-stable processes explicitly char¬ 
acterise spatial dependence, BHM can be more flexible and pragmatic by 
allowing for inclusion of physical mechanisms in terms of covariates and 
random effects. The major issue with BHM is the conditional independence 
assumption of the extremes, whereas for max-stable processes it is model im¬ 
plementation and flexibility. Copula models lie somewhere in between since 
the dependence of the extremes is modelled by the copula assuming that the 
marginal distributions are separable from this dependency structure [Sang 
and Gelfand (2010)]. 

In this paper, we adapt BHM as the modelling framework mainly be¬ 
cause of their flexibility in allowing for (temporal) covariate effects along 
with a versatile spatial dependency structure through random effects. BHM 
generally assume independence of the extremes for given values of the co¬ 
variates and random effects (conditional independence), although they can 
be extended to model spatial extremal dependence by including max-stable 
processes [Reich and Shaby (2012)]. For the application to extra-tropical 
cyclones, we believe conditional independence to be a reasonable working 
assumption. Much of the dependency between successive cyclones has been 
shown to be induced by modulation of rates by time-varying climate modes 
and so can be accounted for by including appropriate covariates [Mailier 
et al. (2006), Vitolo et al. (2009)]. 

There has been recent interest in spatial BHM for extremes since their 
introduction by Casson and Coles (1999). In Cooley, Nychka and Naveau 
(2007) and Cooley and Sain (2010), a CPD and a point process model are 
used to model extreme precipitation where the spatial dependence is charac¬ 
terised by Caussian random effects in the formulation of model parameters 
but without incorporating temporal nonstationarity. Caetan and Crigoletto 
(2007), Heaton et al. (2011) and Sang and Celfand (2009) allowed tempo¬ 
ral structure in BHM through time-varying covariates where the conditional 
model is a CEV distribution. Turkman, Turkman and Pereira (2010) used a 
similar model where the conditional model is a CPD. In this paper, we use 
the computationally efficient MCMC algorithm from Cooley and Sain (2010) 
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based on recent work on Markov random fields [Rue and Held (2005)] and 
add temporal covariates, to account for temporal trends and variations. We 
use the point process model for extremes as the conditional model: it utilises 
more of the data than GEV models and, unlike GPD models, inference is 
invariant to the choice of threshold. 

2.3. Data. Objective feature-identification software [Hodges (1994)] was 
used to extract cyclone tracks from 6-hourly National Genter for Environ¬ 
mental Prediction Glimate Forecast System (NGEP-GFS) re-analysis data 
[Saha et al. (2010)] available over the period 1979-2009. Individual cyclone 
tracks are identified by tracking local maxima in relative vorticity just above 
the boundary layer (about 1.5 km above sea level). The minimum sea-level 
pressure (MSLP) and its location are recorded every 6 hours throughout the 
lifecycle of each cyclone. We use sea-level pressure as a measure of cyclone 
intensity mainly because this variable is well observed and has smooth vari¬ 
ation during the lifetime of a cyclone, unlike other possible variables such as 
wind speed or vorticity. Figure 1(a) shows a map of cyclone tracks defined 
by 6-hourly MSLP recordings for a period with high cyclone activity. Only 
a subset of tracks is plotted: ones where any 6-hourly MSLP value reached 
below 960 hPa. Typical damaging cyclones over Europe reach values in the 
range 940-970 hPa, for instance, Anatol: 953 hPa [Ulbrich et al. (2001)] and 
Kyrill: 962 hPa [Mitchell-Wallace and Mitchell (2007)], whereas the lowest 
ever recorded Braer cyclone reached 913 hPa off the North-West of Scotland 
in January 1993 [Odell et al. (2013)]. 

Although wind speed could also have been used, exploratory analysis 
suggests that extreme MSLP and maximum wind speed are strongly depen¬ 
dent, as to be expected from simple balance arguments. Above the surface 
boundary layer outside equatorial regions, centrifugal and Goriolis forces are 
approximately balanced by the pressure gradient force. Hence, wind speeds 
above the boundary layer in extra-tropical cyclones are proportional to pres¬ 
sure gradients (gradient wind balance). Surface pressure gradients in turn 
are strongly related to the cyclone MSLP since extra-tropical cyclones have 
similar synoptic spatial dimensions (the so-called Rossby scale). Hence, from 
such simple dynamical meteorology arguments, MSLP and maximum wind 
speeds are expected to be extremally dependent and so will convey similar 
information. Let random variables W and Z denote maximum (6-hourly) 
wind speeds at about 1.5 km above the surface (on the 925 hPa pressure 
surface) and negated MSLP (obtained by multiplying MSLP by —1), respec¬ 
tively, with associated 6-hourly recorded values wt and zt- Figure 2(a) shows 
a plot of Wt against zt‘. there is strong positive association with the loess 
smoother indicating a nearly linear relationship. To better visualise extremal 
dependence. Figure 2(b) shows the empirical copula obtained by producing 
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(a) 




Eig. 1. (a) Cyclone tracks for the October 1989 to March 1990 extended winter. Only a 

subset of tracks is plotted: ones where any 6-hourly MSLP value reached below 960 hPa. 
Nadir positions are denoted with solid circles, (b) Sea-level pressure versus latitude and 
(c) latitude for two of the cyclone tracks in (a). 


a scatter plot of the empirical probabilities ql = (rank( 2 ;i) — l)/(n — 1) and 
[Stephenson et ah (2008)], where n is the total number of 6-hourly 
recorded values. This transforms out the margins to uniform distributions 

iz) (tu) 

since ql and ql are estimates of the cumulative distribution functions 
(CDFs) Fz{Z) and Fw{Z). Strong dependence of the extremes is evident 
from the convergence of the points in the upper right-hand corner of the 
graph. 

Figure 2(c) shows estimates of the extremal dependence measure x [Coles, 
Heffernan and Tawn (1999)] defined as x = x(p)) where x{p) = 

Fic{Fz{Z) > p\Fw{W) > p). As p —>■ 1, x{p) Oj implying asymptotic in¬ 
dependence, so we also show x{p), another measure of strength of extremal 
dependence, in Figure 2(d). The quantity x = limp_>.i x(p) measures the 
strength of extremal dependence within the class of asymptotic indepen¬ 
dence. Since x{p) remains positive but does not tend to 1, we conclude that 
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Fig. 2. (a) Wind speed against negated sea-level pressure with an associated loess fit (grey 

line). The intersecting lines are the values —960 hPa and 45 m/s for pressure and wind 
speed, respectively, representing the same high empirical quantile for each variable, (b) 
Empirical copula of wind speed and pressure along with the associated quantile lines from 
(a), (c) The associated extremal dependency measure x(p)> o,nd (d) x(p) P- The 95% 
confidence intervals in (c) and (d) are based on the Normal approximation to proportions 
and are calculated as introduced in Coles, Heffernan and Tawn (1999). 


there is a positive nonasymptotic association at extremes of negated MSLP 
and maximum wind speed, so either variable could potentially be used to 
investigate extremes (see Appendix A.2 for details on x and x). 

Figure 1(b) and (c) show plots of MSLP against latitude and longitude, 
respectively, for two particular cyclone tracks in the 1989-1990 winter [Fig¬ 
ure 1(a)]. The plots illustrate not only the tendency of intense cyclones to 
move in a west-to-north direction but also the fact that MSLP decreases 
(cyclone deepening) as the cyclone propagates in space and time, to reach 
a minimum (which we assume approximates the unobserved value of the 
cyclone nadir) before it starts increasing again until the end of the life cycle. 
Understanding the limiting strength of the nadirs is an important aspect in 
the study of extra-tropical cyclones. However, the rate of growth of cyclones 
depends on the large-scale atmospheric environment that they pass from, 
so the pressure limit of cyclone nadirs will vary with the spatial location of 
the cyclone. By only considering the nadir from each track, we focus on a 
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fundamental limiting property of cyclones, namely how deep they can get 
in general rather than how deep they can get in specific spatial locations. In 
other words, we are interested in spatial variation in cyclone intensity rather 
than maximum local cyclone impact. 

The analysis of nadirs only, also helps to eliminate dependency between 
successive 6-hourly MSLP measures and reduces the amount of data from 
313,557 6-hourly measurements to 17,230 nadirs. Figure 3(a) shows the (re¬ 
analysis) nadir from each track in the Atlantic region where dots in black are 
nadirs with sea-level pressure lower than 960 hPa. However, a single value 
for the threshold defining the extremes is not appropriate and the definition 
of extremeness should vary spatially. For example, a damaging cyclone in 
the Mediterranean is likely to be considered a weak one over Scandinavia. 

Note that we nse cyclone tracks from a reanalysis data set mainly because 
generally cyclone track observations for the extra-tropics are not readily 
available. However, reanalysis data are outpnt from climate models with as¬ 
similated historical observational data. There is much smoothing/ 
interpolation of the observational data when creating a reanalysis data set, 
so the interpretation of any results obtained here is conditional on the effects 
of such smoothing. 

3. Model specification and model fitting. 

3.1. Spatial discretisation. Conventional Bayesian spatial models gener¬ 
ally rely on the assnmption that data are either gridded or they come from 
fixed locations in space [see Banerjee, Carlin and Gelfand (2004)], where one 
or more observations are available at each location. Extreme nadirs, however, 
behave like a spatial marked point process where both location of occnrrence 
and magnitude are random. To utilise such Bayesian models, we propose for 
simplicity to discretise space by imposing a finite grid and to consider the 
minimum possible size A for each grid cell, to ensure that enough data are 
available for estimation in each cell. Inference should not be sensitive to the 
choice of grid spacing, provided it is fine enough (in the limit A —)■ 0 one 
should obtain the original marked point process). Sensitivity analysis for A 
is an important part of the concept (see Section 3.6). 

For spatial marked point processes, estimation is only possible after mak¬ 
ing assumptions about spatial (and temporal) structure. The assumption 
made by discretising is that, conditional on a cell-specific random effect and 
possible covariates, the extreme events (nadirs) within each cell come from 
the same distribntion. The cell-specific random effects are spatially depen¬ 
dent to allow for correlation between events in neighbonring cells. We also 
assnme that the events can occur anywhere within the cells, with equal prob¬ 
ability. Importantly, redefining space into discrete grid cells also provides a 
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Fig. 3. (a) Map of all cyclone nadirs: dots in black represent nadirs deeper than 960 hPa, 

(b) map of recorded X{s,t) that are greater than the threshold (90th empirical quantile) 
in each grid cell and (c) map of thresholds in each cell. 
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way of defining extremes in space: as values below a cell-varying threshold 
or as the r = 1,2,... largest values, in fixed time periods. 

More generally, conditional on a given spatial or spatio-temporal depen¬ 
dence structure between cells, nadirs are modelled using an appropriate ex¬ 
treme value model. This is a hierarchical model where at the top of the 
hierarchy, random effects and covariates define a spatio-temporal process 
which modulates the process, giving rise to extreme nadirs. 

3.2. Spatial grid. Conventionally, extreme value modelling is applied to 
the upper tails so the nadirs are negated to obtain variable X{s,t), where 
s refers to the grid cell and t refers to time. We may think of X{s,t) as 
the depth of a cyclone so that high values of X{s,t) correspond to low 
values of MSLP. We divided the domain in Figure 3(a) into N = 150 5° x 
10° grid cells. The threshold u{s) in each cell was defined as the empirical 
90th quantile of X{s,t). We performed exploratory threshold analysis using 
mean residual life plots [Coles (2001)], ensuring that the 90th empirical 
quantile was an appropriate threshold choice. Figure 3(b) shows the map of 
the extremes (1736 nadirs) and Figure 3(c) shows the map of u{s). Note that 
in Figure 3(c), three cells are highlighted: cells containing coordinates (5.2°E, 
60.2°N), (0°E, 5°N) and (3.5°W, 40.2°N) marked in white crosses. These 
coordinates relate to the cities of Bergen, London and Madrid, respectively, 
and will be used throughout the paper for illustration of results, as they 
adequately span Europe in terms of latitude. 

3.3. Model specification. To model the depth X of negated nadirs, we 
consider the point process model for extremes [Coles (2001)] —conditional 
on spatial random effects and temporal covariates. Eor some high threshold 
u of X, this model is parametrised in terms of the location, scale and shape 
parameters of the GEV distribution, namely, p, a and ^ (see Appendix A.l). 
We use the notation X ~ PF{n,a,^,u). Introducing spatial and temporal 
variation, let X{s,t) be the depth in grid cell s G S at time t G T, where S 
and T are the space and time domains, each a fixed subset of 2-dimensional 
and 1-dimensional Euclidean space, respectively. Extending the approach of 
Cooley and Sain (2010), we model the X{s,t) in the following way: 

(1) X{s,t)\e^{s),fi 2 {s)-^PPin{s,t),a{s,t),({s),u{s)), 

(2) fj.{s,t) = (3^ + /3^zi{t) + /32{s)z2{t) + 9^{s), 

(3) log{a{s,t)) = 13^ + (3^ zi{t) + (s), 

(4) ^(,) = /l« + 0«(,) 

for ip = p,,a,^. Defining vectors = (C/’^(l),..., [/'^(N))' for 

U(s) = lu>^{s),U^{s),U^{s)y and U = (U^^,U^, U?)', the spatial level of 


10 


T. ECONOMOU, D. B. STEPHENSON AND C. A. T. FERRO 


the model is as follows: 

(5) (0^(s),r(s),0«(s))'|U(5)~N(U(5),diag(T)-i), 

(6) U = (U^,U^,U«)'~N(0,n~^), 

(7) /32(s) ~N(z^,(/)2), 

where zi is the latitude of the occurrence and Z 2 is the North Atlantic 
Oscillation (NAO) value (see Section 3.4 about covariate selection). The 
spatial random effects 9^{s), 9'^{s) and 9^s) dehne spatial variability in /r, 
log(cr) and ^ across the cells, after allowing for covariates. The r-year return 
level, that is, the (1 — l/r)th quantile of X(s,t) in cell s and time t, is given 
by 

(8) Xi_i/^{s,t) = fj.{s,t) + ^^^((-log(l - 1/r))"^^*^ - 1). 

As in Cooley and Sain (2010), vectors Ijb are modelled jointly using a 
separable formulation [Banerjee, Carlin and Gelfand (2004), Chapter 7], so 
that the precision matrix is = T 0 W. The matrix T is an unknown 3x3 
positive definite symmetric matrix and W is an x A^ proximity matrix 
dehning spatial proximity between the N cells. Therefore, the dimension of 
ft is 3N X 3N. Here, spatial proximity is based on nearest neighbours so 
that off-diagonal elements of W are wtj = — 1 if cells i and j are adjacent 
and Wij = 0 otherwise, whereas diagonal elements Wi^i = — Wij [see 
Bailey and Gatrell (1995), pages 261-262 for other examples of proximity 
measures]. 

Each vector U^, is modelled by an Intrinsic AutoRegressive (lAR) 

spatial model [Banerjee, Carlin and Gelfand (2004)]. The lAR model uses 
the proximity matrix and a single unknown parameter to control the spatial 
dependency structure (see Appendix A.3). Here, there are three such param¬ 
eters for each of U^, U'^, and they are found in the diagonal of T. (Note 
that the value of r is conventionally fixed beforehand to avoid nonidenti- 
fiability between r and the diagonal of T [Banerjee, Carlin and Gelfand 
(2004)].) Dependence between U^,U°',U^ is modelled using 3 parameters, 
the off-diagonals of T, each controlling the strength of dependence. Allowing 
explicitly for this dependence can aid the MCMG estimation discussed in 
Section 3.5, in terms of convergence to the posterior and also mixing of the 
MCMC samples. 

The NAO parameter /S 2 (s) is spatially variable but in an unstructured 
way. This ensures that /32(s) share information to aid estimation in cells 
with few events but less so compared to using a structured (lAR) spatial 
prior. Parameter reflects the overall NAO effect on /i(s,t). 

We complete the model specihcation by defining the prior distributions of 
the hyperparameters. The intercepts /3 q,/3q,/3q were given Gaussian priors 



MODELLING EXTREME STORMS 


11 


with large variance, and means (/2, logd, ^), calculated as means of indepen¬ 
dent maximum likelihood fits of point process models in each cell. For param¬ 
eters we assumed a flat Gaussian prior with zero mean and large 

variance. The prior distribution 7r(-) for is chosen so that oc 1/^^ 

[Gelman et al. (2013), Chapter 3], whereas for T and P we use a Wishart 
prior with 3 degrees of freedom (uninformative) and a mean that relates to 
the variability of /x, cr and ^ across cells (see Section 3.5). 

3.4. Covariate selection. This was performed by adding explanatory vari¬ 
ables to a “null” model: the model in (l)-(4) without zi or Z 2 . Models were 
compared using the Deviance Information Criterion (DIG), a model selection 
criterion for Bayesian models [Spiegelhalter et al. (2002)] and by investigat¬ 
ing whether posterior distributions of associated parameters are centred at 
zero with relatively large variance. 

The model in (l)-(4) was first implemented with the addition of latitude, 
longitude, latitude squared, longitude squared and an interaction term be¬ 
tween longitude and latitude as covariates in both fJ.{s,t) and log(cr(s, t)). 
This allows for large-scale spatial trends, leaving the local spatial depen¬ 
dence to the random effects. It also relaxes the assumption of complete spa¬ 
tial randomness of extreme events within a cell, both in terms of occurrence 
and intensity. In principle, nonparametric surfaces can also be considered 
for smoothing large-scale spatial trends [see Davison, Padoan and Ribatet 
(2012), page 173 for references], but this was not deemed necessary here. To 
quantify the effect of large-scale climate patterns, two climate indices were 
also considered as covariates: the North Atlantic Oscillation (NAO) and the 
East Atlantic Pattern (EAP), both of which have been shown to be influen¬ 
tial for extra-tropical cyclones [Mailier et al. (2006), Seierstad, Stephenson 
and Kvamsto (2006), Pinto et al. (2009), Nissen et al. (2010)]. No covariates 
were considered for the shape parameter ^(s) since this is a particularly dif¬ 
ficult parameter to estimate, however, it was allowed to vary between cells. 
Out of all possible covariate combinations, the lowest DIG value occurred 
for the particular model formulation in (l)-(4). The posterior distributions 
of “insignificant” parameters (e.g., ones relating to longitude) had means 
and medians very close to zero. 

It is well known that the NAO has influence on the development of extra- 
tropical cyclones [Pinto et al. (2009)]. By definition, the NAO index is stan¬ 
dardised to have zero mean and unit variance, and here it was defined as 
5-day nonoverlapping averages from 1979-2009. Figure 4(a) shows the time 
series of NAO and Figure 4(b) shows the histogram of NAO where the values 
of 2 and —2 are marked, as we consider these as high and low NAO threshold 
values throughout the rest of this paper. Figure 4(c) and (d) show extreme 
values of X{s,t) for which NAO > 2 and NAO < —2, respectively. There is 
a clear North-South pattern in the Central Atlantic, implying NAO has a 
notable effect on extreme cyclones. 
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(c) 
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Fig. 4. (a) Time series of NAO defined as a 5-day average of daily NAO, (b) histogram 

of NAO along with vertical lines marking the values —2 and 2, (c) occurrences of recorded 
nadirs where the associated NAO value was greater than 2 and (d) less than —2. 


3.5. Estimation by Markov chain Monte Carlo. For all ran¬ 

dom effects 9^{s) and /32{s), and parameters (3^ and /3f, were sampled by 
Metropolis-Hastings, specifically using a random walk sampler. The inter¬ 
cepts Pq were sampled from their full conditionals using Gibbs sampling, 
by treating them as intercepts in the mean for each 9'^{s). Samples of 
andT were drawn using Gibbs sampling, utilising the specific 
techniques in Cooley and Sain (2010). Both v and (fP were sampled from 
their full conditionals: Gaussian and scaled inverse-y^, respectively. 

Note that when the lAR model is used as a prior it is improper: the density 
does not integrate to 1. So, to make the intercept terms /3g identifiable, the 
rows of W must sum to zero. This in turn imposes the restriction that 

The parameter t was set to (0.1,10,100)h These values were chosen by 
fitting independent point process models in each cell and investigating the 
level of variability between cells for fi{s), log((T(s)) and C(s), not only to 
reflect the difference in scale for the three parameters but also to make sure 
that most of the variability is modelled by the random effects U^,U'^,U^ 
and not r. If values in r are too small, then the variability in each 9'^{s) is 
forcibly large and may cause problems in estimating the diagonal of T which 
relates to the variability of each Sensitivity analysis was performed to 
ensure these values have little effect on inference (not shown for conciseness). 
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0 2000 4000 6000 6000 10000 0 2000 4000 6000 8000 10000 

MCMCsamplM MCMCnmptM 

Eig. 5. (a) Deviance samples from each of the three MCMC chains. Vertical lines denote 

the burn-in and the total number of simulations. Samples between the two lines are used 
for inference, (b) Samples of the shape parameter ^(s) for the grid cell containing London. 


The Wishart prior for the precision matrix T was given the following 
mean: diag(0.02,4,40)'. As with r, these values were calibrated by fitting 
independent point process models and were chosen to reflect the associated 
levels of variability for each of /u(s), log(cr(s)) and ^(s). 

The model in (l)-(4) was implemented in R [R Development Core Team 
(2012)] using three parallel MCMC chains. These were run on a workstation 
with a 3.07 GHz i7 processor and the processing speed for each chain was 30 
seconds for 1000 samples. A total of 50,000 samples were collected per chain 
and thinned by 5 to reduce auto-correlation. After thinning, the first 3000 
samples from each chain were discarded based on a trace plot of deviance 
(minus twice the log-likelihood) shown in Figure 5(a). Convergence in the 
deviance is a good indication of convergence to the joint posterior of all 
parameters [Gelman et al. (2013)]. Summarising, 21,000 posterior samples 
were used to calculate posterior distribution statistics for the parameters. 
Figure 5(b) shows an example trace plot of .^(s) for the grid cell containing 
the London coordinate. 

3.6. Sensitivity to grid cell size. A purely spatial model [i.e., model (1)- 
(4) without zi and Z 2 ] and a stationary model [i.e., model (l)-(4) without 
zi, Z 2 and the random effects] were implemented for different grid configu¬ 
rations. For each model, the 100-year return level (i.e., the level exceeded 
by the annual maximum in any particular year with probability 0.01) of 
X{s,t) was calculated using (8), for each of the three coordinates marked 
in Figure 3(c). Figure 6 shows the posterior mean of the 100-year return 
level against the number of cells in each grid configuration along with 95% 
credible intervals for each model. Convergence of the return value, as the 
number of cells increases, is evident for the spatial model (although this 
varies slightly due sampling variation). The random effects pool information 
spatially, whereas the stationary model ignores neighbouring cells, resulting 
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Fig. 6. Dots are posterior means of the 100-year return level of X{s,t) versus number 
of cells in different grid specifications, along with 95% credible intervals. Left /(a) and 
(d)/, middle /(b) and (e)/ and right /(c) and (f)/ panels refer to the Bergen, London and 
Madrid cells, respectively. Top /(a), (b), (c)/ and bottom /(d), (e), if)] panels refer to the 
stationary and the spatial models, respectively. For reference, the deepest recorded value of 
X{s,t) in each cell is shown with a cross symbol. 


in failure to converge, especially over London. Pooling also results in notably 
smaller credible intervals for the spatial model—note that the intervals are 
skewed. We chose = 150 cells for the analysis so that all cells have an 
adequate number of nadirs (ranging from 14 to 376). 

4. Results. Posterior distributions for global parameters are summarised 
in Table 1. Latitude has a positive linear effect on both the location and log- 
scale parameters of extreme cyclone depth X(s,t). The overall NAO effect 
u is positive, in agreement with findings from previous studies [Pinto et al. 
(2009)]. To assess MCMC convergence, the Gelman and Rubin R multi¬ 
chain diagnostic was used for each of our model parameters [Gelman et al. 
(2013)]. The R values for each parameter in Table 1 are all close to unity, 
suggesting convergence. 

Figure 7 shows posterior means and standard deviations of /i(s,t), cr{s,t) 
and ^{s). Much of the spatial structure in the extreme nadirs comes from the 
location and scale parameters. The posterior means for the shape parame¬ 
ter are more uniform and generally negative, apart from one cell over 
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Table 1 

Summary of parameter posterior distributions 


Parameter 

Prior 

Posterior mean (s.e.) 

95% Cr.I. 

R 

l3f (Latitude) 

N(0,100) 

4.71 (0.62) 

[3.61,5.94] 

1.13 

Pf (Latitude) 

N(0,100) 

0.12 (0.07) 

[0.00,0.25] 

1.03 

Overall NAO effect v 

N(0,100) 

1.21 (0.24) 

[0.77,1.66] 

1.01 

Variance NAO effect cfp 

OC 

5.6 (1.85) 

[3.09,10.25] 

1.00 


N(-944.1,100) 

-987.4 (0.51) 

[-988.5,-986.5] 

1.04 

PS 

N(5.7,100) 

2.03 (0.06) 

[1.91,2.15] 

1.05 

Pi 

N(-0.19,100) 

-0.13 (0.03) 

[-0.18,-0.07] 

1.10 


Iceland. Exploring this further, the two deepest nadirs in the reanalysis oc¬ 
curred in this cell, and they are considerably lower than the rest of the nadirs 
in the vicinity. A return level plot from the particular cell indicated that the 
two nadirs (one of them being from the record-breaking Braer cyclone) un¬ 
duly influenced the sign of the shape parameter. This has been quantified 
by removing those two points and refitting the model, however, this being 
an analysis of extremes, it makes little sense to remove such values. 

A negative shape parameter implies that the distribution of extreme cy¬ 
clone depth X{s,t), at time t and cell s, has an upper bound given by 
a{s,t)/^{s) — Here this corresponds to a lower limit on nadir sea- 

level pressure. Many of the posteriors for ^(s) do have some mass over the 
positive real line [see, e.g.. Figure 5(b)]. However, except for the Iceland cell, 
the negative masses for ^(s) are all greater than 0.5, therefore, we can use 
the negative posterior ^(s) samples to obtain a conditional posterior distri¬ 
bution for the estimated lower limit. The posterior means of these limits are 
shown in Figure 8(c) for NAO = 0. The limit for the cell containing Bergen is 
890.0 hPa [193.0, 932.6] and for the London cell it is 943.0 hPa [714.8, 959.4], 
whereas in the Madrid cell it is 953.5 hPa [537.9, 978.7]. The 95% credible 
intervals are skewed and noticeably wide, which is to be expected given we 
are trying to estimate the 100th percentile. The lower bounds on some of 
these intervals are too low to be physically plausible and this reflects the 
fact that the statistical model is not constrained by physical mechanisms. 
Note also that there is considerable literature focusing on the problem of 
estimating upper/lower bounds of distributions. See de Haan and Ferreira 
[(2006), Chapter 4] for a detailed discussion and a description of both max¬ 
imum likelihood and moment estimators for bounds arising from extreme 
value distributions. In addition, Einmahl and Magnus (2008) provide re¬ 
fined estimators for bounds of world records in athletics and their respective 
sampling distributions. 

The posterior means and standard deviations of the NAO effects /32{s) are 
shown in Figure 7(g) and (h), respectively. A positive effect is prominent in 
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Fig. 7. Posterior means for (a) ^{s,t), (c) a{s,t), (e) ^(s) and (g) /32(s) and standard 
errors m (b), (d), (f) and (h), respectively, where zi{t) is latitude at centre of grid cell 
and Z 2 {t) = 0. 
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Eig. 8. Estimated lower limits of nadir sea-level pressure for (a) NAO = 2, (b) 
NAO = —2 and (c) NAO = 0. (d) shows the difference between (a) and (b). 


the area where cyclones deepen the most: in the vicinity of Iceland, northern 
Europe and Scandinavia. A negative effect is also apparent, effectively over 
Spain and the Azores. This North-South NAO effect in the central Atlantic 
is consistent with the exploratory diagnostics in Figure 4(c) and (d). Maps 
of the estimated lower limit for NAO = —2 and NAO = 2 are given in Fig¬ 
ure 8(a) and (b). To better see the effect of NAO on the estimated lower 
limit. Figure 8(d) shows the difference in hPa between the estimated lower 
limits for NAO = 2 and NAO = —2. The difference can get up to 25 hPa in 
the area where NAO has the biggest effect, that is, northern Europe and 
Scandinavia. 

Figure 9 shows return level plots of X(s, t) for the Bergen-London-Madrid 
grid cells, for NAO = ±2. Note that this is not a goodness-of-fit test, as 
each point in these plots (the recorded value) is associated with a different 
NAO value, whereas the return level curves are calculated at NAO = ±2. 
A positive/negative NAO effect is noticeable in the Bergen/Madrid cells, 
confirming the NAO North-South effect. No NAO effect is evident in the 
London cell. The horizontal line in each plot is the estimated cyclone depth 
limit, suggesting that for all three cells, nadirs could have been much deeper 
than the ones recorded. 
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Fig. 9. Individual grid cell return level plots (posterior means) with 95% credible in¬ 
tervals. Observed values shown in solid circles. Top panel: NAO = 2; bottom panel: 
NAO = —2. Left panel: Bergen cell; middle panel: London cell; right panel: Madrid cell. 
Horizontal lines are estimated upper bounds of X{s,t) for NAO = 2 (top) and NAO = —2 
(bottom). 

Therefore, we also consider the quantity vr(s,f) = Pr(X(s,t) > Xm{s)), 
where Xm{s) is the negated minimum recorded nadir in grid cell s for the 
30-year period. (Note that this is equivalent to describing how unusual the 
recorded depth was, rather than the probability of ever getting deeper than 
the recorded 30-year minimum nadir.) We transform the GEV parameters to 
reflect the distribution of 30-year, rather than yearly depth values: a = aS^ 
and jl = iJL + d(l — where S = 30. Figure 10(a) shows 7r(s, t) for values 

of NAO associated with Xm{s). There are high values of 7r(s,f), especially 
over western Europe. Figure 10(b) shows 7r{s,t) for NAO = 2, indicating 
that for a positive NAO phase there is high probability of deeper nadirs over 
Europe, Iceland and Scandinavia. For NAO = —2, Figure 10(c) attributes 
high probability of deeper nadirs over Spain, Portugal, west of France and 
also over the Azores region. Furthermore, Figure 10(d) shows the difference 
in hPa between the estimated depth limit for MSLP [Figure 10(a)] and Xm{s) 
for each cell. For most cells, the difference is in the range of 10-50 hPa, 
while for cells over Iceland the range is 80-110 hPa, indicating the 30-year 
reanalysis is not long enough to capture nadir depths near the estimated 
limits. 

We use posterior predictive checking Gelman et al. [(2013), Ghapter 6] to 
assess model fit. This compares each observation, x(s,f)obs) to the poste¬ 
rior predictive distribution for replications, Ai(s,f)rep, of X{s,t) given the 
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Fig. 10. Probability of observing a deeper nadir than the recorded 30-year deepest nadir 
in each cell: (a) Calculated for NAO values associated with the recorded values nadirs, (b) 
NAO = 2 and (c) NAO = —2. (d) The difference in hPa between the estimated depth limit 
and the deepest recorded 30-year nadirs in each cell. 


data, D, used to fit the model. If the observations do not behave as if they 
are sampled from their posterior predictive distributions, then this indicates 
poor model fit. Samples of (s, t)rep were obtained by simulating from GEV 
distributions with parameters equal to draws from their joint posterior distri¬ 
bution and then the posterior predictive means and 95% posterior predictive 
intervals were approximated from these samples. We plot the observations 
of (a) the deepest 30-year nadirs and (b) the deepest yearly nadirs against 
the corresponding posterior predictive means and intervals in Figure 11(a) 
and (b), respectively. None of the observations seem extreme with respect 
to the posterior predictive distributions: the 45-degree line falls well within 
the prediction intervals. 

We also calculate the probability integral transform (PIT), z{s,t) = 
Pr(X(s,t)rep < a^('S,t)obs \ D), of each observation relative to its posterior 
predictive distribution. If the model is a good fit, then the z{.s,t) should 
follow a uniform distribution on the interval (0,1). For each grid cell, s, we 
plot the probability points {i — l)/{n{s) — 1) for i = 1,... ,n{s) against the 
order statistics of the z{s,t) values for that cell, where n(s) is the number 
of observations in cell s. Departures from the 45-degree line indicate poor 
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Fig. 11. Recorded versus predicted values of: (a) 30-year deepest nadirs in each cell and 
(b) yearly deepest nadirs in each cell. The predicted values are the means of the posterior 
predictive distributions while the grey shaded area represents the associated 95% prediction 
intervals. 


model fit. We indicate the sampling variation that would be expected in 
these plots when the model is perfect by pointwise 95% confidence intervals, 
constructed by simulating samples of size n{s) from the uniform distribu¬ 
tion on (0,1). Figure 12 shows these plots for Bergen, London and Madrid. 
No points fall outside the 95% intervals, indicating adequate fit. Note that 
PIT values are often used in forecast verification; see, for instance, Gneiting, 
Balabdaoui and Raftery (2007) and references therein. Although histograms 
are the more conventional way of displaying PIT values, here we only have 
a few data points for each cell, so we use probability-probability plots. 

5. Conclusions. We have implemented a flexible model, adapted from 
Cooley and Sain (2010), to reanalysis cyclone data in what we believe to 
be the first study that simultaneously models both the spatial and tempo¬ 
ral structure of extreme extra-tropical cyclones. Using (1) spatial random 



PIT value 


PIT value 


PIT value 


Fig. 12. Probability-probability plots of theoretical Unif(0,1) probabilities versus proba¬ 
bility integral transform (PIT) values z{s,t) for the Bergen, London and Madrid cells. The 
95% confidence intervals reflect sampling uncertainty. 
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effects, (2) latitude as a covariate and (3) a 150 cell spatial regularisation, 
spatial variation was adequately modelled in the extremal behaviour of the 
cyclones. The North Atlantic Oscillation was nsed as a covariate and was 
found to have a significant effect on extremal cyclone behaviour, especially 
over Northern Europe and the Iberian peninsula. 

Although this is a hrst step toward studying the spatio-temporal be- 
havionr of extreme cyclones, the analysis relies on assumptions which may 
oversimplify the problem: (1) the creation of an artificial grid, (2) the choice 
of threshold in each cell and (3) the subjective choice of spatial proximity. 
The choice of the grid is a potential weakness which can introdnce bias, as 
both the nnmber of cells and their shape are subjectively chosen. Techniques 
such as Dirichlet tessellation or Delaunay triangulation [Illian et al. (2008)] 
may be usefnl for defining a more optimal “data-driven” grid. The shape of 
the cells is particularly important if one is interested in modelling data along 
cyclone tracks rather than individual points as in onr application. If interest 
was in the relative spatial cyclone impact, one could use cell-specific rather 
than cyclone-specific nadirs, rendering the rectangular cells inappropriate. 
Hexagonal cells wonld be more appropriate as illnstrated in an application 
to tropical cyclones in Eisner, Hodges and dagger (2012). Threshold choice 
in each cell may also prove to be an issue. Ideally, model fit should be one 
of the criteria for choosing the threshold. Eor the application in this paper, 
three different thresholds were considered: the 85%, 90% and 95% quantile 
in each cell. Model fit diagnostics (Eigures 11 and 12) indicated worse fit for 
the 85% quantile, and an identical fit for the higher quantiles—which is why 
we selected the 90% quantile for model implementation. To avoid choosing 
the threshold altogether, one might instead estimate the threshold from the 
data. Eor example, we have explored the possibility of using a mixtnre model 
as in Erigessi, Haug and Rue (2002) where the threshold is estimated but 
which also allows all available data (not jnst the extremes) to be used for 
each cell, which in tnrn allows the use of a finer grid. Last, the proximity 
structure used to define the covariance matrix of the spatial random effects 
is also an assnmption which can affect the degree of spatial smoothing. The 
spatially random occnrrence of cyclone nadirs was “marginalised” here by 
dividing the region into grid cells, whereas one should ideally try to model 
both the spatial occurrence and intensity at the same time, for example, 
by using spatial marked point process models. Nevertheless, despite these 
assnmptions, the model in this study is flexible enough to be used in other 
similar studies, for example, ones involving tropical cyclone wind speed max¬ 
ima or cyclone-related peak precipitation. 

APPENDIX 

A.l. Point process model for extremes. The point process model for 
extremes involves a bivariate variable Y = {X,T), with T G [0,1] being a 
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scaled random variable associated with time and X G M a random variable 
associated with intensity. The model is a marked point process, which for 
X > u [a. high threshold) under some linear normalisation and mixing cri¬ 
teria [Smith (1989)] behaves like a nonhomogeneous Poisson process with 
intensity function 


(9) 


X{x,t) 


1 

a 


1+C 


X — /i 
a 


- 1 / 5-1 


provided that 1 -|- (^/cj)(x — /i) > 0. The intensity function X{x,t) is zero for 
I + {C/cr){x — fi) <0. The exceedance rate is explicitly modelled in terms of 
the mean number of exceedances in the time interval [ti,t 2 ]- 


A([ti,t 2 ] X {u,oo)) = {t 2 -h) 


1 + e 


n — /t 
a 


-i/S 


The likelihood given observations yi = {xi,ti) in region [0,1] x (u, cx)) is 


( 10 ) 

( 11 ) 


L(/i,(T,^; x,t) = exp|—riy y J A(x,t)dxdt| JjA(xj,tj) 


= exp< —Uy 


1 + ? 


u — 
a 


-i/« 




where Uy is the number of years of observed data so that parameters /r, 
a and ^ correspond to the GEV distribution of yearly maxima. Because 
the time variable T does not actually appear in (9) and thus in (11), we 
use the concise notation X ~ PP{p,a,^,u) as in Section 3.3. The likelihood 
contribution from a single event {xi,ti) is 

L(/x, a, u) = exp<^ [t* - U-i] 


1 + ^ 


u — y 
a 




X{xi,ti) 


for i = 0,..., n where n is the number of events. Note that to = 0 and that 
the likelihood contribution, for the time interval between the last event oc¬ 
currence and t = 1, is the probability of no events in the interval, that is, 


exp 


^ ^?/[l 


1 + 1 


u — y, 
a 


-i/« 


The conditional model in (l)-(4) was implemented using the likelihood 
(11) for each cell. However, because of the temporal covariates, the outermost 
integral over time in (10) is impossible to calculate analytically unless one 
knows explicitly how the covariates evolve in time. A remedy is to approx- 
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imate the integral: divide the time range in small intervals with endpoints 
0 = /ci, A: 2 , ..., A:j = 1 and assume the function is constant in each interval. 
The integral 

dt 

is thus approximated by the Riemann sum 


/ 


l + Cis) 


u{s) - 
a{s,t) 


1 

J 


J r 


E 

i=l 


i+as) 


u{s) - 
a{s,ki) 


-i/m 


where J is the number of intervals. In practice, J is determined by observa¬ 
tions of the covariates for all data (not just the extremes). 


A.2. Measures of extremal dependence. The measure of extremal de¬ 
pendence 0 < X < 1 between random variables Z and W is defined as 

X= lini Fr{Fz{Z) >p\Fw{W) > p) = 

p—>■! 

where Fz and Fw are the respective distribution functions of Z and W. 
The other extremal dependence measure — 1 < X E 1 is defined as 

X=lim-21ogPr(Fz(^)>;,)-i = lin.j(p). 

^ P^i logPr(Fz(^) > p,Fwiw) > p) 

If X > 0 and x = 1; the two variables are asymptotically dependent and x 
measures the strength of that dependence. If x = 0 and x < 1; the two vari¬ 
ables are asymptotically independent, in which case x measures the strength 
of dependence—within the class of asymptotically independent variables. 
Roughly, X measures the “speed” at which x(p) approaches zero. Coles, 
Heffernan and Tawn (1999) advocate the use of both x and x as indicators 
of extremal dependence, providing complementary information on different 
aspects of that dependence. 


A. 3 . Intrinsic AutoRegressive priors. Consider a grid with N cells. If the 
random effect (j) = (</>(!),..., (j){N)y is assumed to have an lAR prior, then 
0 ~ A^(0, (tW)“^), where W is the adjacency matrix and the conditional 
distribution for each <p{s) given the rest is given by 

4>{s)\(t){-s) ~ iv(0(s), —, 

\ Tm[s) J 

where is 0 excluding (()(s); ^{s) is the average of (t>{—s) that are 

adjacent to (p{s) and m{s) is the number of those adjacencies. 
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